
# Add Eastern HDD variables ----------------------------------------------------
# Load NARR degree-day data
narr_dt = file.path(dir_rds, 'degreeDayNARR.rds') %>% readRDS()
# Convert date variable to date format
narr_dt[, date := ymd(date)]
# Key datasets
setkey(sub_dt, hh_id, begin_date)
setkey(narr_dt, date)
# Merge datasets
sub_dt %<>% merge(
  x = .,
  y = narr_dt[region == 'east', c('date', str_subset(names(narr_dt), '^tot_')), with = FALSE],
  by.x = 'begin_date',
  by.y = 'date',
  all.x = TRUE,
  all.y = FALSE
)
# Clean up
rm(narr_dt); invisible(gc())
# Create five HDD bins (based upon quintiles)
hdd_q = sub_dt[, quantile(tot_hdd_lag_month08, probs = seq(0, 1, 0.2))]
q1 = hdd_q[2]
q2 = hdd_q[3]
q3 = hdd_q[4]
q4 = hdd_q[5]
sub_dt[, `:=`(
  tot_hdd_lag_month08_bin1 = fcase(
    between(tot_hdd_lag_month08, 0, q1), tot_hdd_lag_month08,
    tot_hdd_lag_month08 > q1, q1
  ),
  tot_hdd_lag_month08_bin2 = fcase(
    tot_hdd_lag_month08 < q1, 0,
    between(tot_hdd_lag_month08, q1, q2), tot_hdd_lag_month08 - q1,
    tot_hdd_lag_month08 > q2, q2 - q1
  ),
  tot_hdd_lag_month08_bin3 = fcase(
    tot_hdd_lag_month08 < q2, 0,
    between(tot_hdd_lag_month08, q2, q3), tot_hdd_lag_month08 - q2,
    tot_hdd_lag_month08 > q3, q3 - q2
  ),
  tot_hdd_lag_month08_bin4 = fcase(
    tot_hdd_lag_month08 < q3, 0,
    between(tot_hdd_lag_month08, q3, q4), tot_hdd_lag_month08 - q3,
    tot_hdd_lag_month08 > q4, q4 - q3
  ),
  tot_hdd_lag_month08_bin5 = fcase(
    tot_hdd_lag_month08 < q4, 0,
    between(tot_hdd_lag_month08, q4, Inf), tot_hdd_lag_month08 - q4
  )
)]

# Estimates with Eastern HDD instrument ----------------------------------------
# Marginal price: Polynomial, 2nd degree
hdd_poly2 = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_mrg_lag2) ~ 
      poly(tot_hdd_lag_month08, 2, raw = TRUE) +
      poly(tot_hdd_lag_month08, 2, raw = TRUE):socal,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Marginal price: Polynomial, 3rd degree
hdd_poly3 = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_mrg_lag2) ~ 
      poly(tot_hdd_lag_month08, 3, raw = TRUE) +
      poly(tot_hdd_lag_month08, 3, raw = TRUE):socal,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Marginal price: Polynomial, 4th degree
hdd_poly4 = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_mrg_lag2) ~ 
      poly(tot_hdd_lag_month08, 4, raw = TRUE) +
      poly(tot_hdd_lag_month08, 4, raw = TRUE):socal,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Marginal price: Polynomial, 5th degree
hdd_poly5 = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_mrg_lag2) ~ 
      poly(tot_hdd_lag_month08, 5, raw = TRUE) +
      poly(tot_hdd_lag_month08, 5, raw = TRUE):socal,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Marginal price: Polynomial, 6th degree
hdd_poly6 = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_mrg_lag2) ~ 
      poly(tot_hdd_lag_month08, 6, raw = TRUE) +
      poly(tot_hdd_lag_month08, 6, raw = TRUE):socal,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Marginal price: 5 bins
hdd_bins = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_mrg_lag2) ~ 
      tot_hdd_lag_month08_bin1 + tot_hdd_lag_month08_bin1:socal +
      tot_hdd_lag_month08_bin2 + tot_hdd_lag_month08_bin2:socal +
      tot_hdd_lag_month08_bin3 + tot_hdd_lag_month08_bin3:socal +
      tot_hdd_lag_month08_bin4 + tot_hdd_lag_month08_bin4:socal +
      tot_hdd_lag_month08_bin5 + tot_hdd_lag_month08_bin5:socal,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Table of results
etable(
  hdd_poly2, hdd_poly3, hdd_poly4, hdd_poly5, hdd_poly6, hdd_bins,
  style.tex = style.tex('aer'),
  fitstat = ~ r2 + n, tex = TRUE,
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'pooled-mrg-hdd-city.tex'),
  replace = TRUE
)



